In the theory of differential equations concepts related to stability are very well established and understood. Parallel theories for numerical schemes in Euclidean spaces are also well-developed. From the beginning, the notion of absolute stability, using a linear test equation was the prevalent tool. Later, Dahlquist, Burrage and Butcher were leading the development of a numerical stability theory that also makes sense for nonlinear problems, that of B-stable or non-expansive methods, set in Hilbert spaces. In more recent times, it has become popular to construct and analyze numerical schemes for differentiable manifolds, entirely intrinsic methods, examples are the Lie group integrators. For the study of nonlinear stability, Riemannian manifolds, or even Finsler manifolds, seem to be a useful framework for studying non-expansive methods. Building on the work by Kunzinger et al. (2006) and Simpson-Porco and Bullo (2014), we shall suggest a way to generalize the notion of B-stability to Riemannian manifolds. We study some simple integrators designed for simulating dynamical systems on manifolds, in particular, we study methods based on geodesics. We show some unexpected results for the geodesic version of the implicit Euler method where it turns out that its nonlinear stability depends on the sectional curvature of the metric. We shall also discuss the existence and uniqueness of solutions to the nonlinear system that must be solved in each time step for implicit methods. Explicit schemes can not be unconditionally non-expansive, but we shall see that under a special type of monotonicity condition they can be non-expansive when the step size is appropriately bounded. The theory we present has important applications for damped mechanical systems evolving on Lie groups and homogeneous spaces, such as in multibody dynamics. But recently, it has also found applications in dynamical systems-based neural networks. When the data input to the neural networks is manifold-valued, it becomes an issue to design architectures that have the 1-Lipschitz property in order to improve robustness.